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Abstract 



The experiment KASCADE observes simultaneously the electron-photon, muon, 
and hadron components of high-energy extensive air showers (EAS). The analysis of 
EAS observablcs for an estimate of energy and mass of the primary particle invokes 
extensive Monte Carlo simulations of the EAS development for preparing reference 
patterns. The present studies utilize the air shower simulation code CORSIKA with 
the hadronic interaction models VENUS, QGSJet and Sibyll, including simulations 
of the detector response and efficiency. By applying non-parametric techniques the 
measured data have been analyzed in an event-by-event mode and the mass and 
energy of the EAS inducing particles are reconstructed. Special emphasis is given to 
methodical limitations and the dependence of the results on the hadronic interaction 
model used. The results obtained from KASCADE data reproduce the knee in the 
primary spectrum, but reveal a strong model dependence. Owing to the systematic 
uncertainties introduced by the hadronic interaction models no strong change of 
chemical composition can be claimed in the energy range around the knee. 

Key words: cosmic rays; energy spectrum; mass composition; knee; EAS 
PACS: 96.40.De 
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1 Introduction 



The basic astrophysical questions in high-energy cosmic rays (CR) relate to 
the sources, the acceleration mechanisms and the propagation of CR through 
space. In particular, the observation of the change of the power law slope 
(the knee [1]) of the all-particle spectrum at an energy of a few times 10^^ cV 
has induced considerable interest and experimental activities. Nevertheless, 
despite of many conjectures and attempts, the origin of the knee phenomenon 
has not yet been convincingly explained. 

Due to the rapidly falling intensity and low fluxes, cosmic rays of energies 
above 10^'^ eV can be studied only indirectly by observations of extensive air 
showers (EAS) which are produced by successive interactions of the cosmic 
particles with nuclei of the Earth's atmosphere. EAS develop in the atmo- 
sphere as avalanche processes in three different main components: the most 
numerous electromagnetic (electron-photon) component, the muon component 
and the hadronic component. The properties of EAS are usually measured 
with large ground-based detector arrays. In most experiments only one or two 
components are studied. The KASCADE experiment [2,3] studies all three 
main components simultaneously and a large number of shower parameters 
are registered for each event. Their analysis to determine the properties of 
the primary particle are obscured by the considerable fluctuations of EAS 
development. 

The analysis of the EAS variables to deduce the properties of the primary 
particle relies on the comparison with Monte Carlo simulations (MC) of the 
shower development (sec Figure 1), including the detector response. Usually 
only one or two EAS parameters are measured and various simplified proce- 
dures are used to describe the relation between the observed EAS properties 
and the nature and energy of the primary particle. The simplification often im- 
plies the use of parameterizations of the average behavior, which may bias the 
results and limit the accuracy because fiuctuations are neglected or not prop- 
erly accounted for. For the analysis of multivariate parameter distributions 
and accounting for fluctuations more sophisticated methods are needed. The 
decades-old Bayesian methods and the neural network approaches, currently 
in vogue, meet these necessities. The methods facilitate an event-by-event 
analysis. 

In the present paper we report on an investigation of the energy spectrum and 
mass composition of cosmic rays in the energy range of 10^^ — 10^^ eV, based on 
the analysis of 700,000 EAS events. A subset of approximately 8000 showers 
with cores near the center of the hadron calorimeter yields information on all 
three components and has been studied in more detail. Following the analysis 
scheme shown in Figure 1, the simulated showers calculated with the simula- 
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Fig. 1. Twofold way of the EAS analysis procedure. 

tion program CORSIKA [4] have been convoluted with the apparatus response 
using the GEANT code [5]. Non-parametric procedures [6] yield not only an 
estimate of the primary energy and mass composition, but they also allow to 
specify the uncertainty of the results in a quantitative way. In addition, we 
specify the dependence of the results on the hadronic interaction models. The 
necessity to invoke such models in an energy range beyond the experimental 
limits of accelerator experiments, implies a model dependence of the results 
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on the energy spectrum and mass composition. Quantifying this model de- 
pendence is one of the objectives of the present paper. The model dependence 
is illustrated by using two different interaction models for the analysis. The 
dependence imphes not only the degree to which a particular EAS observable 
is correlated to energy and mass of the primary particle, it shows also how 
sensitively different EAS observables reveal primary mass. As an example, the 
mass composition depends on the particular set of observables being consid- 
ered simultaneously in the analysis if the model is inconsistent with the data 
in all internal correlations. 

It should be stressed that the present study emphasizes the methodical aspects 
of how to infer energy spectrum and mass composition of CR rather than 
providing a final answer. This would require improved statistical accuracy 
both in experiment and simulation and, first of all, a reduction of systematic 
uncertainties due to the incomplete knowledge of high energy interactions. 
Nevertheless our findings on spectrum and mass composition are compatible 
within the methodical accuracy to the results of other experiments. 



2 The KASCADE experiment 



The detector installation of the experiment KASCADE (KArlsruhe Shower 
Core and Array DEtector) [2,3] is located on the site of the Forschungszcntrum 
Karlsruhe, Germany (8° E, 49° N; 110m a.s.l.). The three major components 
of the detector system (Figure 2) are 

• an Array of scintillation detectors, 

• a Central Detector: an arrangement of several different detector compo- 
nents, basically a hadron iron sampling calorimeter using liquid ionization 
chambers and 

• a Muon Tracking Detector (MTD) using limited streamer tubes. 

The Array covers an area of about 200 x 200 m^ and consists of 252 detector 
stations. These are organized in 16 clusters and placed on a square grid of 
13 m separation. The detector stations contain liquid scintillation counters 
(e/7 detectors) of 0.79 m^ area each and plastic scintillators of 0.81 m^ each 
(yU detectors; E^^'^'^^ = 230 MeV), the latter covered by a shielding of 10 cm lead 
and 4 cm steel. The inner four clusters (60 stations) contain four e/7 detectors 
per station but no /i detectors while the outer 12 clusters (192 stations) have 
two e/7 detectors and four /i detectors per station. The reconstruction of 
the EAS data measured with the Array provides the basic information about 
lateral distributions and total intensities of the electron-photon (shower size 
Nq) and muon components {N^; see Section 4), the location of the EAS core 
and the direction of incidence. 
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Fig. 2. Schematic layout of the KASCADE experiment. 

The layout and performance of the Central Detector are described in refer- 
ence [7]. The finely segmented hadron calorimeter is the main part of in the 
Central Detector system. It consists of a 20 x 16 m^ stack of about 4000 tons of 
iron with eight horizontal gaps. The calorimeter thickness corresponds to 11 
interaction lengths Aj for vertical hadrons. The detectors, measuring the en- 
ergy deposit of the traversing charged particles, are ionization chambers filled 
with the room temperature liquids tetramethyl-silane (TMS) or tetramethyl- 
pentane (TMP), inserted into the gaps of the iron stack and read out by 40,000 
electronic channels. From their signals the impact point, the direction and the 
energies of individual hadrons are reconstructed. In particular, the number of 
EAS hadrons with energies larger than 100 GeV (A^jf ^^°'"~^°^), the energy of 
the most energetic hadron observed in the shower {E^'^^) and the energy sum 
of all reconstructed hadrons (J2 Eh) are deduced as shower observables (see 
Section 4). 
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Table 1 

KASCADE detector components used in the present analysis. 



Detector 


Total Area 


Threshold E^n 


Observables 


array e/7 


490 m^ 


5MeV 




array fj, 


622 m^ 


230 MeV X sec^ 




MWPCs 


122 m^ 


2.4 GeV X sec 6' 


at;, D_6, Dg 


calorimeter 


320 m^ 


50GeV 





A layer of 456 scintillation detectors, each with a size of 0.45 m^, is mounted 
in the third gap at a depth of 2.2 Aj. It is used for triggering the Central 
Detector system, for muon detection (with a threshold of EJ^^'^^ — 490 MeV), 
and to determine arrival time distributions [8] . 

In the basement of the central building, below the iron stack and 77 cm of 
concrete, two layers of multi-wire proportional chambers (MWPCs) are ar- 
ranged as a tracking hodoscope, covering an area of 122 m^ [9] . The MWPCs 
register muons with an energy threshold Ej^^^^ — 2.4 GeV and provide the ob- 
servable A*"*, i.e. number of reconstructed muons in the MWPCs. Due to the 
good position resolution, the MPWCs register also the spatial distribution of 
the high-energy muons together with traversing secondaries produced in the 
absorber by high-energetic hadrons, whose pattern has been shown to carry 
valuable information about the mass of the primary particle [10], expressed 
by particular parameters {D_q, Dq) in terms of a fractal moment analysis (see 
Section 4). Information on specific detector details is compiled in Table 1. 



3 Simulations 

The simulations of the EAS development, along the requirements of the analy- 
sis scheme of Figure 1, have used the air shower simulation program CORSIKA 
(vers. 5.62) [4]. The code incorporates several options of high-energy interac- 
tion models and is continuously under improvement. In particular, we consider 
the latest versions of VENUS [11], QGSJet [12] and Sibyll (vers. 1.6) [13]. 
VENUS and QGSJet are models based on the Gribov-Regge theory, and ex- 
trapolate the interaction features in a well defined way into energy regions 
which are far beyond energies available by accelerators, and especially into the 
extreme forward direction. Sibyll is a minijet model used as a hadronic inter- 
action generator in the MOCCA [14] and the AIRES codes [15]. We use it here 
only for demonstration purposes. For the low-energy interactions CORSIKA 
includes the GHEISHA code [16]. The infiuence of the Earth magnetic field 
on charged particle propagation is taken into account. As density profile of 
the atmosphere the U.S. standard atmosphere is chosen [4]. 



7 



Samples of at least 2000 proton and iron-induced showers have been simulated 
with all three models. Additionally for VENUS and QGSJet the intermediate 
mass primaries He, O and Si have been simulated. The energy distribution 
follows a weighted power law with a spectral index of —2.7 in the energy 
range of 10^^ eV to 3.16-10^^ eV, calculated in eight intervals. The zenith angles 
are distributed in the range [13°, 22°]. The centers of the showers are spread 
uniformly over an area which exceeds the surface of the hadron calorimeter by 
2 m on each side. In addition, roughly the same number of simulated events 
with the centers of the showers within the Array are used. The signals observed 
in individual detectors are determined by tracking all secondary particles down 
to observation level and passing them through a detector response simulation 
program based on the GEANT package [5]. 



4 Event reconstruction and selection 

The reconstruction of the EAS observables which is described in detail in 
preceding publications of the KASCADE collaboration [10,17-20], apphes an 
iterative procedure for reconstructing the shower size parameters. In a first 
step the shower core location is determined by a center-of-gravity technique 
from the energy deposit signals of all e/7 counters, and the shower direction 
is estimated by a simple plane fit using the timing information of the Array 
detectors. In addition, as rough first approximations, the electron size A^e and 
muon size are estimated from summation of detector signals, taking into 
account the actual shower core position on the grid. These parameter values 
are initial values for the further reconstruction steps. In the second step the 
shower direction is determined by fitting a conical shape of the shower disc to 
the arrival times of the charged particle component, registered with the e/7 
counters. The lateral distributions and their shape parameters are estimated, 
and and A^e are determined. 

The muon size Nj^ is the muon content within a range of distances from 
the shower core between 40 m to 200 m which is the range accessible to the 
KASCADE experiment (the so-called truncated muon number) [20] . The lower 
limit is chosen to exclude contributions of the electromagnetic and hadronic 
punch-through near the center of the showers. The upper limit corresponds 
to the geometrical acceptance of the KASCADE layout. Figure 3 displays the 
variation of A^*'' and Ne versus the primary energy E, as inferred from EAS 
simulations. Due to various serendipitous features of the KASCADE layout 
IgA^'*'' proves to be nearly proportional to IgE and turns out to be almost 
independent of primary mass (Figure 3 left). This is in contrast to the electron 
size Ng, which exhibits a strong mass dependence for fixed Nj^ as shown in 
Figure 4 (right). 
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5 5.5 6 6.5 7 7.5 5 5.5 6 6.5 7 7.5 



Ig (E/GeV) Ig (E/GeV) 

Fig. 3. The mean values of the truncated muon number Njf and electron number 
A^e vs. primary energy as inferred on basis of the indicated interaction models. For 
sake of clarity only QGSJet predictions are fitted by a linear function in Ig-lg scale, 
in order to emphasize the much more pronounced mass dependence of the shower 
size A'e. 

Contributions to the detector signals of other particles than electrons and 
muons are eliminated by applying a lateral energy correction function to ap- 
propriate particle densities, which are fitted with a likelihood function to the 
Nishimura Kamata Greisen (NKG) formula [21,22]. Values of the radius pa- 
rameters of 89 m and 420 m for electrons and muons, respectively, are used [17]. 
For showers whose cores are located within 91m from the Array centerQ, the 
reconstruction uncertainty is about 2 m for the location of the shower center, 
0.5° for the angle of incidence, and less than 10% and 20% for and 
values, respectively, at primary energies larger than 10^^ eV. 

Muon tracks observed with the MWPCs, reconstructed from pairs of hits in 
the two MWPC layers (vertically separated by 38 cm [10]), are summed up to 
obtain N*. A limit for the reconstructed angle of ±15° in zenith and ±45° in 
azimuth with respect to the shower axis determined from the Array is imposed 
(the azimuth cut is not applied for showers with zenith angles of < 10°). The 
analysis of the number and spatial distribution of the muons and of produced 
secondaries in terms of two generalized multi-fractal dimensions D_q and Dq 
is discussed in [10]. These parameters characterize the spatial distribution 
of muons and high-energy (punch-through) hadrons as well as the degree of 
fluctuations of particles in the shower core. 

The reconstruction of the hadronic shower variables applies appropriate pat- 
tern recognition algorithms [18,19]. Energy clusters found in different detector 

^ This number 91 m results from the extension and the grid spacing of the detector 
array. 
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Fig. 4. The variation of the mean values of the electron shower size (left) and 
number of reconstructed hadrons ]\[^>^00GeV (^j-^gj^t) with A^^"". The predictions of 
different models for proton and iron induced showers are compared with results of 
the measurements. 

layers are traced from lower layers to the topmost one to form a particle track. 
Additionally, the angle of incidence of the track can be deduced by the same 
procedure starting at lower layers, patterns of cascades have to form clusters 
from the remaining energy bunching up to showers according to the already de- 
termined direction. Furthermore, signals in the first layer from the top are not 
used for energy determination, because their electromagnetic punch-through 
distorts the hadron signals. The signals, weighted by the overlying absorber 
thickness are summed up and converted to hadron energies [7]. Similar to the 
shower size Nc, the reconstructed number of hadrons (jV^>iooGeV'j exhibits a 
strong mass dependence for fixed as well (Figure 4; right). 

For the investigation of the primary energy spectrum and mass composition, as 
well as of particular correlations of observables, two sets of data are compiled. 
One set - further referred to as selection I - uses the information of electrons 
and muons from the Array stations only. It allows analysis the data with 
good statistical accuracy, but includes only little information provided by the 
Central Detector. A second data set - henceforth referred to as selection II - 
includes observables measured in the Central Detector, but this data sample 
comprises only a small amount of registered showers. Selection I comprises 
720,000 EAS events, accumulated in 226 days, with primary energies larger 
than E fa and with a maximum core distance from the center of the 

Array of 91 m and with angles of incidence in the range of [13°, 22°] . Selection II 
comprises approximately 8000 high-energy, central showers selected by cuts on 
Nj^ {IgN^ > 3.5), on the core location (-Rcorc < 5m from the center of the 
Central Detector) , with at least one reconstructed hadron of an energy above 
100 GeV and 10 muons observed in the MWPCs. 
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5 Non-parametric analyses 



The present analysis of mass composition and energy spectrum avoids the 
bias inherent in parametric procedures and is performed for individual events 
by use of multivariate non-parametric Bayesian and neural network decision 
methods. In this way we are able to specify, in a transparent and coherent way 
how conclusive and trustworthy our results are, as expressed by true classifi- 
cation and misclassification matrices of the results. A brief outline and more 
details of the applied methods are given in Appendix A and in reference [23]. 

The combination of the total muon content and the shower size has 
been shown to be sensitive to primary mass and is applied in numerous exper- 
imental studies, using suitable parameterizations of the predicted Ig Nfj_/ Ig A^e 
relation with the primary mass. However, as indicated above, the total muon 
content A^"^, although displaying some dependence on primary mass, is a quan- 
tity not easily accessible experimentally without additional assumptions about 
the shape of the lateral muon density distribution at large distances from the 
shower core. Therefore, we prefer to consider the truncated muon number Nj^, 
which - on average - proves to be nearly independent from primary mass (see 
Figure 3), but it is, on the other hand, a rather sensitive energy identifier. 
Thus, at fixed A^*^, the information about the mass is essentially provided by 
the shower size Ae [20]. In cases of other EAS observables mass and energy 
sensitivities are, in general, less well marked, and in principle, each shower vari- 
able carries information simultaneously on mass and energy in a way which is 
additionally affected by the considerable fluctuations of the shower develop- 
ment. The most sensitive EAS observables, Ae and A'*'', display the smallest 
intrinsic and sampling fluctuations. 

5.1 Mass composition 

Due to the limited number of simulated EAS and the correspondingly limited 
statistical accuracy it is hardly reasonable to use the full set of observables 
simultaneously to achieve a reliable result about mass composition (curse of 
dimensionality condition; see Appendix A). Hence we consider simultaneously 
only a few observables. 

Each simulated or measured event is represented by an observation vector 
X — (Ag) A^'^) • • •) of the n observables. Applying the technique described in 
Appendix A the likelihood (probability density distribution) p{x\cui) of an 
event for each class Ui G {p, O, Fe} can be calculated, i.e. the probability of 
an event x belonging to a given class cui. 

As an example, the superposition of the estimated probability density distri- 
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Fig. 5. Superposition of three probability density distributions Ylf=iPi^\'^i)/'^ de- 
duced from QGSJet simulations using the observables N,. and N^^ (left) as well as 
A'^ and A'^^ (right). Events in the dark shaded area mark the region classified as 
iron, middle grey as oxygen, and light grey as proton {selection IT). 

butions, referring to two sets of different observables, are displayed in Fig- 
ure 5 (based on QGSJet simulations). The regions where p(a;|a;p), p(a;|ci;o) and 
p^xlui-pe) are larger than the other two possibilities are colored light, middle 
and dark grey, respectively. The left graph shows the density distribution cal- 
culated in the two dimensional space of the observables and A^'^''. A rough 
separation can be recognized, but also a strong overlapping of the likelihood 
distributions has to be admitted. 

The right-hand graph of Figure 5 shows an example of two observables {N* 

and A''*''), which exhibit only weak mass-discrimination power. Correspond- 
ingly, the density distributions of the three particle types are intermixed, and 
reliable conclusions could not be drawn. In case of selection //the mass compo- 
sition is reconstructed for different sets of observables using the Bayes theorem 
(Equation A.l). When the estimated posterior probability 'p{uji\x) is larger 
than p(ujj\x), then the event is assigned to class cjj, otherwise to class ujj. 
Taking into account (by Equation A.l) the estimated number of incorrectly 
classified events (i.e. misclassification rates) (Table 2) the true proportions of 
the different particle types are reconstructed. 

The classificat (see Appendix: Equation A.l and Fig- 

ure A.l) give the fraction of correctly, Pjj, and wrongly, Pj^, classified events 
with i 7^ j, an example for three mass class is given in Table 2. Of course, 
the sum of each row has to be 100%. In the most probable cases the different 
particle types are identified correctly, but the knowledge of the incorrectly 
classified events could be used for a correction due to the mis-classification. 
In addition, the rates for the intermediate mass particle types. He and Si, 
are given. Helium is mostly classified as protons (57%) and silicon as oxygen 
(54%). Due to the stronger fluctuations and weaker correlations with mass 



12 



Tabic 2 

Classification matrices for three classes (p, O and Fe) and two different models. In 
addition to the classification rates of p, O and Fe, the rates of classified intermediate 
groups He and Si respectively, are given. The used observables are Njf and Ne 
(3.6 < IgiV*"^ < 3.9). 

QGSJc^t VENUS 



Pujj^uji [%] =P =0 uji =Fe cji =p Ui =0 LOi =Fe 





=P 


77 ±3 


21 ±3 


2± 1 


78 ±3 


21 ±2 




LVj 


=He 


57 ±3 


39 ±3 


4± 1 


64 ±3 


32 ±2 


4± 1 


Uj 


=0 


14 ±2 


61 ±3 


25 ±3 


15 ±2 


61 ±4 


24 ±3 


Uj 


=Si 


3±1 


54 ±3 


43 ±2 


3±2 


51 ±3 


46 ±2 


UJj 


=Fe 


1±1 


17 ±2 


82 ±3 




20 ±3 


80 ±3 



Fe Fe 

P((OFel'<)='' P(0>Fel'<)='' 




Fig. 6. Estimated posterior probabilities p{u)i\x) of measured events deduced from 

QGSJet simulations using the observation vector x = {Nq,N^) (left) as well as 
x= (iV;J,Are>iooGeV^^^^) ^^.^-^^-^ 

and/or energy, different sets of observables result in lower true-classification 
rates Pa. In general, if the rates Pa are less than 50%, it is no more possible to 
deconvolute the true proportions by matrix inversion of Pij (Equation A. 5), 
since the matrix P^j becomes singular, signaling that the determination of a 
class uJi is just haphazard. Therefore it is not meaningful to consider more 
than 3 classes, since this would require an analysis of further observables si- 
multaneously, with a number of Monte Carlo simulations larger than presently 
available. 

As a cross-check the estimated postenor probabilities p{uji\x) of a given mea- 
sured event x, belonging to class Ui, can be calculated (Figure 6). The center 
of the triangles shown correspond to equal probability of belonging to any 
class, reflecting the fact that it is nearly impossible to classify the measured 
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Fig. 7. The relative abundances of the classes p, O and Fe vs. IgNj^'^, reconstructed 
on the basis of two different hadronic interaction models and using the EAS observ- 
ables A^e and Nj^. The right graph shows the corresponding mean logarithmic mass 
{In A) vs. IgN^. Statistical (thick lines) and methodical (thin lines) uncertainties 
are indicated as error bars. 



event, while points in the corners satisfy the relation p{uji\x) = 1, i.e. the cor- 
responding event belongs to class cUj with probability unity. Hence, from given 
measured events we obtain information about the probabilities belonging to 
class LOi just as well. Evidently, the set A^e and allows to determine a well 
defined mass composition. In contrast, the set comprising 
J2 Eh is less suitable for mass discrimination mainly due to strong fluctuations 
of the observables as can be seen in Figure 6, right side. 

The results of the composition determination, using the observables and 
A'*'', are given in Figure 7 as relative abundances versus A^'' and mean log- 
arithmic mass (In A) vs. A'*''. The mean logarithmic mass {In A) cannot be 
calculated unambiguously because only abundances of groups of elements can 
be determined. We have assigned (Ap+He) = 2.5, {Aq) = 16, (v4si+Fe) = 42 
to the p, O and Fe group, respectively. This procedure is of course to some 
extent arbitrary, but this is always implicit, when (ln74) is used. The statisti- 
cal errors (thick lines) are calculated according to a multinomial distribution. 
The thin error bars correspond to a methodical uncertainty, calculated by the 
bootstrap method (see Appendix A). It reflects the influence of the limited 
number of Monte Carlo events. Apparently the use of the VENUS model re- 
sults in a lighter composition, compared with the use of QGSJet. A tendency 
towards a heavier composition above the knee (IgA'*'' = 4.15 corresponding 
to 4PeV) is indicated, albeit the statistical and systematic uncertainties do 
not allow a definite conclusion and the results are clearly compatible with an 
energy independent composition. 
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Fig. 8. Mean logarithmic mass {In A) resulting from the analysis of different sets of 
observables vs. IgA'^'^ (QGSJet prediction). The sets displayed on the right do not 
include the observable A'e- The error bars are omitted to simplify the presentation 
of the synopsis. The statistical errors are of the same order of magnitude as given 
in Figure 7. But the systematic errors are larger by 40% — 50% for the data on the 
right graph due to a weaker correlation with mass. 

Results of several other combinations of observables by using the QGSJet 
model are summarized in Figure 8. In general, the tendencies are the same. 
Remarkably, all sets omitting the electron size A'e (right graph) result in a 
heavier composition and a more pronounced increase above the knee. As the 
electron size has the strongest mass sensitivity, as well as the smallest fluctu- 
ations, the mass compositions are predominantly determined by iVg and Nj^ 
(left). Compositions resulting from sets of less sensitive observables differ from 
these values (right). The tendencies are quite similar for the VENUS model, 
but the absolute values are shifted towards a lighter composition as expected 
from Figure 7. 

The fact that different combinations of observables taken into account in the 
analysis, lead apparently to different mass compositions (shown in Figure 8), 
reveals inadequacies of the reference model, i.e. that the degree of the in- 
trinsic correlations of different observables differs from those of the real data. 
Otherwise the determined mass compositions should be identical within the 
statistical errrors. 



5.2 Energy spectra 



To estimate the primary energy E the most important parameters are 
and Nj^ again, where now Nj^ carries most of the information. As data basis 
we use selection I. Due to the large computing time requirements we do not 
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Fig. 9. The relative error vs. primary energy for different classes: Results obtained 
with a network trained with QGSJet samples (left) and the result of the same 
network analyzing VENUS samples (right). 

apply the Bayesian algorithms here and use instead neural networks only. 
In principal there are no basic arguments to prefer one particular method. 
Previous publications have demonstrated the consistency and equivalence of 
neural network and Bayesian methods in EAS analyzes [24,25]. The neural 
networks employed have typically a simple net topology 2x5x2x1, but several 
other topologies have been used to estimate the methodical error of this special 
choice. For training the network (according to Equation A. 7) two independent 
samples have been generated to allow a validation of the results (Appendix A) . 
Before data analysis, the response and the biases of the trained neural network 
have to be carefully scrutinized. For this, the performance of mass classification 
and energy estimation has been investigated with the other MC sample as 
input. We consider the relative deviation of the reconstructed energy E^st from 
the true value E'true, which is known for the simulated samples, more precisely, 
the distribution of AE/E = {E^st — -E'true)/-E'true) whose mean value and the 
standard deviation represent the bias and the energy resolution (relative error) 
of the reconstruction, respectively. 



Figure 9 displays the relative error of the estimated energy for different pri- 
mary particles. The relative error of a network, trained with QGSJet samples, 
is shown in the left part. In general, the bias for the various classes is less than 
3-5%, but the energy resolution (spread) proves to be strongly mass depen- 
dent. As expected, the iron class has the smallest energy spread (ape ~ 21% 
versus ap ~ 38%). The network trained with VENUS samples leads to the 
same results. Also shown (right) are the results of a network, which has been 
trained with QGSJet samples analyzing events generated with the VENUS 
model. With increasing energy the QGSJet trained network underestimates 
systematically the true energy of the VENUS samples. Moreover, the bias 
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Fig. 10. Differential energy spectra resulting from the analysis of data of the 
K AS CADE experiment using two differently trained networks (by VENUS and 
QGSJet samples). The reconstructed energies are compared on an event-by-event 
basis and their differences are given in the inset as relative error vs. the energy 
reconstructed on the basis of the QGSJet model. 



appears to be mass dependent, implying that the degree of the correlation 
between energy and the analyzed shower observables is varying differently for 
the different primary particles and models. This is a caveat to the estimate of 
true energies of the measured event samples, namly that hidden mass depen- 
dent correlations lead to an all-particle spectrum depending on the true mass 
composition. 

Figure 10 presents the reconstructed energy spectra of measured data result- 
ing from the analysis using two different networks, trained with QGSJet and 
VENUS samples. Apparently, the VENUS trained network results in a steeper 
spectrum as compared with the QGSJet findings. It should be emphasized 
that the network used takes into account not only the absolute values of the 
observables Ne and A'"*'', but also their correlations. In order to specify the 
relative error arising from the model dependence, mean value and spread of 
{AE/E)^^^^^ = (^VENUS - ^QGSJet) /^QGSJet are additionally given (inset). 
The variation of this model error displays a change at higher energies, which 
might indicate a change of the composition. 
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Table 3 

Comparison of parameters of the energy spectrum, derived on basis of on the 
VENUS and QGSJet simulations, respectively. The first error is the statistical one. 
The second error represents the systematic uncertainty resulting from the small 
number of simulated event. 

QGSJet VENUS 

71 2.77±0.003±0.03 2.87±0.003±0.04 

72 3.11±0.02 ±0.06 3.25±0.02 ±0.06 
£;knee [10^ GeV] 5.5 ±0.2 ±0.8 4.5 ±0.3 ±0.9 
X^/dof 0.95 1.94 
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Fig. 11. The differential all-particle cosmic ray flux obtained here compared to the 
results reported by Tibet AS7 [29], Akeno [30], and CASA-MIA [27]. The data 
points are multiplied by {E / GeV)^'^ . Only statistical errors are presented. 

The resulting spectra are fitted by the trial function (see Figure 10) 

m-io-(^]"^ ■[i + i^] 1 (1) 



knee 




which accounts for a smoothly changing power law spectrum [26]. The param- 
eter e controls the width of the transition region, and the knee position -Ekncc 
is defined by the center point of the transition region. Asymptotically I{E) 
approaches to power law functions £'~')'1'T2. The parameter uncertainties have 
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been studied by calculating the errors I{E) ± AI(E) using the sampling cor- 
relation matrix. But the resulting error bands are so narrow that it docs not 
visibly differ from the I{E)-lme. The best-fit results are given in Table 3, in- 
cluding statistical errors as well as the methodical error derived from different 
training parameters of the neural network. It is obvious that the statistical er- 
rors are considerably smaller than the systematic uncertainties resulting from 
the small number of simulated events and from interaction models. 

Figure 11 compares the spectra of Figure 10 with results reported by other 
experiments. All measurements, independently from each other, show a steep- 
ening above a particular energy: the knee. But the absolute intensity of the 
flux and the position of the knee obviously differ. This is most likely due to 
different model assumptions or energy conversion functions used. The consid- 
erable deviation between CASA-MIA [27] and the other two experiments may 
be explained in this way. CASA-MIA used the Sibyll model for constructing 
an energy estimator E = f{Nc,Nfj) from the electron and muon sizes. The 
fact that Sybill predicts significantly lower values of and larger values for 
Ne [28] as compared to QGSJet and VENUS, could lead to a systematic shift 
of the spectrum towards lower energies. In view of the considerable model 
dependence of our results, the overlap with some of the other experiments 
should not be taken as evidence for or against any of them. 



5.3 Combined analysis of energy and mass 

As an example. Figure 12 shows the relative abundances of a mass classification 
into three categories, as well as the corresponding mean logarithmic mass, 
resulting from the analysis of central showers (selection II) vs. the estimated 
energy. The observables A^e and Nj^ are used as input parameters. Again, 
the VENUS model leads to a hghter mean logarithmic mass in the considered 
energy range as compared to the QGSJet model. Besides the Bayesian method 
a neural network analysis was performed additionally. The network results are 
denoted by NN. Within the statistical errors the mass composition resulting 
from both pattern recognition procedures agree. Our data reveal a mixed 
composition, becoming lighter when approaching the knee and heavier above 
the knee. This feature appears to be somehow mysterious, but the tendency 
is supported by recent results of the CASA-BLANCA [31] and HEGRA [32] 
experiments. In fact, there are also astrophysical arguments for a minimum of 
the mean mass in the range of the knee [33] . 

In the present status of our analysis procedure it is hardly possible to intro- 
duce more than three classes for the reconstruction of the mass composition. 
If this were to be attempted additional observables had to be included. A finer 
binning of the energy scale (beyond the energy resolution {AE/E)est) for the 
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Fig. 12. Relative abundances reconstructed by Bayes classification vs. the recon- 
structed energy based on the QGSJet model and using A^e and Nj^. Additionally, 
the corresponding mean logarithmic mass (in A) (right; Bayes) and the correspond- 
ing variation resulting from the neural network analysis (NN) arc given. The error 
bars represent the statistical (thick) and methodical (thin) uncertainties. 

spectra of single masses would require to deconvolute the resolution effects. In 
the actual analysis this step has not been performed and only a few represen- 
tative values of the varying mass composition (and no detailed energy spectra 
of the different mass classes) have been presented. To analyze the data beyond 
this limit we need, in the simplest case, to construct from the misclassification 
matrices a matrix ^aa';EE' deconvoluting mass and energy resolution effects. 
Stressing once more the curse of dimensionality (see Appendix A), a very large 
number of simulated events is required for the determination of such a matrix 
(at least 150000 simulated events are needed). For the same reasons we are 
presently unable to infer any significant fine structure from the all-particle 
energy spectrum beyond the resolution {AE /E)est- 

Different sets of observables lead to different mass compositions. This feature 
may indicate that not only the absolute values, but also the correlations of the 
observables are not described satisfactorily by any of the models. Nonetheless, 
a reduced model dependence can be observed when a well marked relation is 
projected out from the correlations of the multivariate distribution. Figure 13 
displays on the left-hand side the relation between electron size A^"c and energy 
E, which shows no difference between simulated and measured samples. Of 
course, this is not surprising, since the pattern recognition tool is just trained 
in that way, such that deviations, incompatible with the statistical accuracy, 
would cast some methodical doubts on the used algorithms. More remarkable 
is the agreement of jVjf ^^^"^^^ vs. primary energy, found by the same network 
though with larger fluctuations of the mean values. That may be explained by 
the reduced mass sensitivity of TVjf and the dominance of the Ne-Nj^ 
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Fig. 13. The projected relations Ig A^e and ig JVjf ^^^'^^^^j respectively, vs. lg(£^/GeV) 
from two neural networks, trained to estimate the energy and mass of the measured 
events using Ne and Nj^ as EAS observables. 

correlation (compare the observable sets in Figure 8 left). Nevertheless, within 
the statistical significance level (in terms of hypotheses tests hke student t- 
test) no difference between data and model predictions can be stated. 



6 Discussion and conclusion 



The present paper aims at presenting methods of a determination of primary 

energy spectrum and mass composition of cosmic rays in the energy range 
lO"*^^ — 5 ■ 10^^ cV by an event-by-cvcnt analysis of EAS data. The specific me- 
thodical feature is the use of a non-parametric approach, studying multivariate 
distributions of a number of EAS observables [34,25]. 

The present approach to obtain information about the EAS primaries has 
following merits: 

• It specifies the inevitable model dependence of any statement about spec- 
trum and mass composition, introduced through the patterns provided by 
the Monte Carlo simulations on basis of a particular hadronic interaction 
model. 

• The model dependence is not only revealed by the results from the analysis 
of single EAS observables when comparing different hadronic interaction 
models, but the approach specifics also the degree of correlations between 
different observables used for the multivariate analysis. 

• This feature provides the possibility to test a specific hadronic interaction 
model by exploring the internal consistency of the results, when the outcome 
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of different sets of observables are considered. This aspect is of greatest 
importance for approaching the best model reproducing the observations in 
the most consistent way. 
• Comparing the KASCADE findings with other experiments shows that the 
discrepancies between results can well be attributed to the different inter- 
action models employed. 
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A Non-pcirametric statistical inference 



Pattern recognition techniques are efficient tools to determine the correct asso- 
ciation of a given sample to a certain category or class. Prom the measurements 
or simulations of a physical phenomenon, a set of quantities (observables) is 
obtained, like or Nc, which defines an observation vector x. This obser- 
vation vector serves as the input to a procedure based on decision rules, by 
which a sample is assigned to one of the given classes. Thus it is assumed that 
an observation vector is a random vector x whose conditional density function 
p{x\u!i) depends on its class cui (e.g. p, O and Pe classes). 

In the following we consider so called non-parametric techniques like Bayes 
classifiers and artificial neural networks [6] . The term non-parametric indicates 
that the representations of the distributions (like probability density functions 
of Bayes classifiers or weights of neural networks) are no more specified by 
a-priori chosen functional forms. They are constructed through the analysis 
process by the given data distributions themselves. 

It should be immediately emphasized that there are some important limita- 
tions. In case of a finite set of random samples, the dimension n of the random 
vector X is limited by the following condition: When considering each compo- 
nent of an n-dimensional observation vector by M divisions, the total number 
of cells is M" and is increasing exponentially with the dimensionality of the 
input space. Since each cell should contain at least one data point this require- 
ment implies that the size of training samples (or reference pattern samples) 
needed to specify the non-parametric mapping, is increasing correspondingly. 
This condition is called the curse of dimensionality [35] and prohibits the 
simultaneous (multivariate) analysis of a larger number of EAS observables, 
when the size of training samples is too small. 



A.l Bayesian decision rule 



The Bayes classifier is a powerful algorithm but time consuming with large 
memory requirements. However, its performance is generally excellent and 
asymptotically Bayes optimal, so that the expected Bayes error (see below) is 
less than or equal to that of any other technique [36]. The estimated probability 
densities converge asymptotically to the true density with increasing sample 
size [37,38]. 

The method is based on the Bayes Theorem [39] 

plxluji) X P(u!i) likelihood x prior 

p{u;i\x) = <^ posterior = — (A.l) 

p{x) normalization factor 
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with p{x) = Z^j=iP(a;|a;j)P(a;j), which holds if the different N hypotheses cui 
(i.e. classes) are mutually exclusive and exhaustive. By a prior and a normal- 
ization factor the theorem connects the likelihood for an event a; of a given 
class uji with the probability of a class cuj, being associated to a given event 
X. The prior gives the a priori knowledge of the relative abundance of each 
class and is major basis of debates on Baycsian inference procedures. It is 
nearly always the best to follow the advice given by Bayes himself [39], gen- 
erally known as Bayes ' Postulate (occasionally also referred to as Principle of 
Equidistribution of Ignorance) : So far there exists no further knowledge, the 
prior probabilities should be assumed to be equal 

1 ^ 

-'^ i=l 



In the fortuitous case that the likelihood functions p{x\uJi) are known for all 
populations, the Bayes optimal decision rule is to classify x into class cui, if 

p{u!i\x) > p{u!j\x) (A. 3) 



for all classes ujj ^ LJi, as illustrated (with the mis-classification probabilities) 
in Figure A.l. 

To construct an estimate p{x\u!i) of the likelihood p{x\ui) of class Wj, the k-th 
simulated event Xki is assumed to have a sphere- of -influence where it con- 
tributes to the probabihties (see Figure A.l). There are various procedures to 
specify these contributions whose superpositions lead to continuous likelihood 
functions, replacing the frequency distributions of discrete simulated events 
Ni of each class uJi in the n-dimensional observation space. A standard choice 
of such spheres are multivariate normal distributions, because they are simple, 
well behaved, easily computed and have been shown in practice to perform 
well: 



p{x\uji) = 




(A.4) 



The Mahalanobis m,etric \ \x — X}^i\ \ = [x — Xki)'^C~^{x — x^i) is used, because 
the observables are transformed to unit variances by the sampling covariance 
matrix Cj for each class uJi, resulting in equal importance of all components. 

The scaling parameter a controls the width of the sphere- of- influence and is 
obtained by the median of an ordered statistics of estimated Bayes errors for 
different initial values of a [23] . The Bayes error e represents the total sum 
(integral) of mis-classified events and is given in case of two classes by the 
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Fig. A.l. Schematic illustration of the construction of two one-dimensional (over- 
lapping) likelihood functions p{x\LOp^Fe}, approximated by Gaussian distributions 
{sphere- of -influence) for each event, indicated on the abscissa (left). Classification 
using the Bayes decision showing the proportion of mis-classified events by the 
hashed areas (right) {P{uJFe) = P{^p)) ■ 

simple relation e = / m\ii{p{uji\x),p{(jj2\x)} ■p{x)dx (hatched areas in Fig- 
ure A.l right). To account for the mis-classification, the rates P, 



i.e. the probability of an event x G cuj being classified in the class uj, are esti- 
mated by the leave-one-out method (also called jack-knifing). Each simulated 
event is held back once while the others are used to estimate the association 
of this particular event. By a, so called, bootstrap method different subsets of 
each simulated class are used to perform the leave-one-out method to give an 
asymptotically unbiased estimate of the variance of the Pij [36]. Thus the true 
number of events n* can be deduced from the classified events nj by a matrix 
inversion: 



ypr 

j 



with = P. 



(A.5) 



A. 2 Neural networks 

An artificial neural network can be considered as a nonlinear transfer function 

f -.W — >W (A.6) 



mapping a bounded euclidian space of dimension p to another space of dimen- 
sion q. The, so called, multilayer feed-forward neural network is organized in 
L different layers: an input layer, L — 2 hidden layers, and an output layer. 
Each layer I consists of a certain number rii of units (neurons), which carry 



27 



on the signals to the next layer. The, so called, network topology specifies the 
number of units in each layer: ni x n2 x . . . x n^^i x n^. An output unit 
of the output layer L is determined for each observation vector {input units) 
and class cui entering the input layer and should be close to the true value 
tki, given by the labeled simulation events in terms of a well defined measure. 
Thus, the error function E{w) 

I N ^ Ni 

= 9 E 7^ E iVmLiXki, W) - t,,)' (A.7) 
^ i=l k=l 

has to be minimized. For each layer /, except of the input, the outcome of 
each neuron m is calculated by a weighted sum of the output of neurons of 
the last preceding layer. Additionally an activation function f{z) is applied to 
the sum 

Vmi = f{z) = / (E <-i • + ■ (A.8) 

A convenient practice is to use the Fermi function f{z) = 1/(1 + exp(— 2)). 
The most common algorithm for the network training, i.e. minimizing E(w) , is 
an adjustment of the weights wf]_i and by a stochastic minimization pro- 
cedure [23] or alternatively by the, so called, back-propagation algorithm [40]. 
There exist different other algorithms or extended versions of this basic back- 
propagation, which try to circumvent problems in finding the global minimum 
or sticking in a local minimum. Additional problems arise, if the training pro- 
cess leads to an overtraining of the network by adopting the properties of the 
training samples, but cannot give satisfactory results, when it is applied to 
another validation set. Thus, in a generalization phase one has to control the 
quality of the network with an independent labeled set of samples. 

In general, the output UmL is a continuous function. Hence not only the classi- 
fication can be done applying neural networks, but also parameter estimation 
(regression) is possible, e.g. the estimation of the primary energy of EAS 
events. In case a classification into N classes is required, the output ymL of 
the network is divided into N regions, each representing a single class [23] . 

In previous publications the consistency and equivalence of neural network and 
Bayes classifier results in EAS analysis have been demonstrated [24,25]. The 
classification rates inferred from both procedures do not differ significantly. 
Thus, an adequate choice of the particular decision rule and of the appropriate 
algorithm is just a matter of the actual conditions like computing time and 
memory workload. 
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